home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / poly / test.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  16.0 KB  |  464 lines
  1. /* poly/test.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <gsl/gsl_math.h>
  23. #include <gsl/gsl_test.h>
  24. #include <gsl/gsl_ieee_utils.h>
  25. #include <gsl/gsl_poly.h>
  26.  
  27. int
  28. main (void)
  29. {
  30.   const double eps = 100.0 * GSL_DBL_EPSILON;
  31.  
  32.   gsl_ieee_env_setup ();
  33.  
  34.   {
  35.     double x, y;
  36.     double c[3] = { 1.0, 0.5, 0.3 };
  37.     x = 0.5;
  38.     y = gsl_poly_eval (c, 3, x);
  39.     gsl_test_rel (y, 1 + 0.5 * x + 0.3 * x * x, eps,
  40.           "gsl_poly_eval({1, 0.5, 0.3}, 0.5)");
  41.   }
  42.  
  43.   {
  44.     double x, y;
  45.     double d[11] = { 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1 };
  46.     x = 1.0;
  47.     y = gsl_poly_eval (d, 11, x);
  48.     gsl_test_rel (y, 1.0, eps,
  49.           "gsl_poly_eval({1,-1, 1, -1, 1, -1, 1, -1, 1, -1, 1}, 1.0)");
  50.  
  51.   }
  52.  
  53.   /* Quadratic */
  54.  
  55.   {
  56.     double x0, x1;
  57.  
  58.     int n = gsl_poly_solve_quadratic (4.0, -20.0, 26.0, &x0, &x1);
  59.  
  60.     gsl_test (n != 0, "gsl_poly_solve_quadratic, no roots, (2x - 5)^2 = -1");
  61.   }
  62.  
  63.   {
  64.     double x0, x1;
  65.  
  66.     int n = gsl_poly_solve_quadratic (4.0, -20.0, 25.0, &x0, &x1);
  67.  
  68.     gsl_test (n != 2, "gsl_poly_solve_quadratic, one root, (2x - 5)^2 = 0");
  69.     gsl_test_rel (x0, 2.5, 1e-9, "x0, (2x - 5)^2 = 0");
  70.     gsl_test_rel (x1, 2.5, 1e-9, "x1, (2x - 5)^2 = 0");
  71.     gsl_test (x0 != x1, "x0 == x1, (2x - 5)^2 = 0");
  72.   }
  73.  
  74.   {
  75.     double x0, x1;
  76.  
  77.     int n = gsl_poly_solve_quadratic (4.0, -20.0, 21.0, &x0, &x1);
  78.  
  79.     gsl_test (n != 2, "gsl_poly_solve_quadratic, two roots, (2x - 5)^2 = 4");
  80.     gsl_test_rel (x0, 1.5, 1e-9, "x0, (2x - 5)^2 = 4");
  81.     gsl_test_rel (x1, 3.5, 1e-9, "x1, (2x - 5)^2 = 4");
  82.   }
  83.  
  84.   {
  85.     double x0, x1;
  86.  
  87.     int n = gsl_poly_solve_quadratic (4.0, 7.0, 0.0, &x0, &x1);
  88.  
  89.     gsl_test (n != 2, "gsl_poly_solve_quadratic, two roots, x(4x + 7) = 0");
  90.     gsl_test_rel (x0, -1.75, 1e-9, "x0, x(4x + 7) = 0");
  91.     gsl_test_rel (x1, 0.0, 1e-9, "x1, x(4x + 7) = 0");
  92.   }
  93.  
  94.   {
  95.     double x0, x1;
  96.  
  97.     int n = gsl_poly_solve_quadratic (5.0, 0.0, -20.0, &x0, &x1);
  98.  
  99.     gsl_test (n != 2,
  100.           "gsl_poly_solve_quadratic, two roots b = 0, 5 x^2 = 20");
  101.     gsl_test_rel (x0, -2.0, 1e-9, "x0, 5 x^2 = 20");
  102.     gsl_test_rel (x1, 2.0, 1e-9, "x1, 5 x^2 = 20");
  103.   }
  104.  
  105.   /* Cubic */
  106.  
  107.   {
  108.     double x0, x1, x2;
  109.  
  110.     int n = gsl_poly_solve_cubic (0.0, 0.0, -27.0, &x0, &x1, &x2);
  111.  
  112.     gsl_test (n != 1, "gsl_poly_solve_cubic, one root, x^3 = 27");
  113.     gsl_test_rel (x0, 3.0, 1e-9, "x0, x^3 = 27");
  114.   }
  115.  
  116.   {
  117.     double x0, x1, x2;
  118.  
  119.     int n = gsl_poly_solve_cubic (-51.0, 867.0, -4913.0, &x0, &x1, &x2);
  120.  
  121.     gsl_test (n != 3, "gsl_poly_solve_cubic, three roots, (x-17)^3=0");
  122.     gsl_test_rel (x0, 17.0, 1e-9, "x0, (x-17)^3=0");
  123.     gsl_test_rel (x1, 17.0, 1e-9, "x1, (x-17)^3=0");
  124.     gsl_test_rel (x2, 17.0, 1e-9, "x2, (x-17)^3=0");
  125.   }
  126.  
  127.   {
  128.     double x0, x1, x2;
  129.  
  130.     int n = gsl_poly_solve_cubic (-57.0, 1071.0, -6647.0, &x0, &x1, &x2);
  131.  
  132.     gsl_test (n != 3,
  133.           "gsl_poly_solve_cubic, three roots, (x-17)(x-17)(x-23)=0");
  134.     gsl_test_rel (x0, 17.0, 1e-9, "x0, (x-17)(x-17)(x-23)=0");
  135.     gsl_test_rel (x1, 17.0, 1e-9, "x1, (x-17)(x-17)(x-23)=0");
  136.     gsl_test_rel (x2, 23.0, 1e-9, "x2, (x-17)(x-17)(x-23)=0");
  137.   }
  138.  
  139.   {
  140.     double x0, x1, x2;
  141.  
  142.     int n = gsl_poly_solve_cubic (-11.0, -493.0, +6647.0, &x0, &x1, &x2);
  143.  
  144.     gsl_test (n != 3,
  145.           "gsl_poly_solve_cubic, three roots, (x+23)(x-17)(x-17)=0");
  146.     gsl_test_rel (x0, -23.0, 1e-9, "x0, (x+23)(x-17)(x-17)=0");
  147.     gsl_test_rel (x1, 17.0, 1e-9, "x1, (x+23)(x-17)(x-17)=0");
  148.     gsl_test_rel (x2, 17.0, 1e-9, "x2, (x+23)(x-17)(x-17)=0");
  149.   }
  150.  
  151.   {
  152.     double x0, x1, x2;
  153.  
  154.     int n = gsl_poly_solve_cubic (-143.0, 5087.0, -50065.0, &x0, &x1, &x2);
  155.  
  156.     gsl_test (n != 3,
  157.           "gsl_poly_solve_cubic, three roots, (x-17)(x-31)(x-95)=0");
  158.     gsl_test_rel (x0, 17.0, 1e-9, "x0, (x-17)(x-31)(x-95)=0");
  159.     gsl_test_rel (x1, 31.0, 1e-9, "x1, (x-17)(x-31)(x-95)=0");
  160.     gsl_test_rel (x2, 95.0, 1e-9, "x2, (x-17)(x-31)(x-95)=0");
  161.   }
  162.  
  163.   {
  164.     double x0, x1, x2;
  165.  
  166.     int n = gsl_poly_solve_cubic (-109.0, 803.0, 50065.0, &x0, &x1, &x2);
  167.  
  168.     gsl_test (n != 3,
  169.           "gsl_poly_solve_cubic, three roots, (x+17)(x-31)(x-95)=0");
  170.     gsl_test_rel (x0, -17.0, 1e-9, "x0, (x+17)(x-31)(x-95)=0");
  171.     gsl_test_rel (x1, 31.0, 1e-9, "x1, (x+17)(x-31)(x-95)=0");
  172.     gsl_test_rel (x2, 95.0, 1e-9, "x2, (x+17)(x-31)(x-95)=0");
  173.   }
  174.  
  175.   /* Quadratic with complex roots */
  176.  
  177.   {
  178.     gsl_complex z0, z1;
  179.  
  180.     int n = gsl_poly_complex_solve_quadratic (4.0, -20.0, 26.0, &z0, &z1);
  181.  
  182.     gsl_test (n != 2,
  183.           "gsl_poly_complex_solve_quadratic, 2 roots (2x - 5)^2 = -1");
  184.     gsl_test_rel (GSL_REAL (z0), 2.5, 1e-9, "z0.real, (2x - 5)^2 = -1");
  185.     gsl_test_rel (GSL_IMAG (z0), -0.5, 1e-9, "z0.imag, (2x - 5)^2 = -1");
  186.  
  187.     gsl_test_rel (GSL_REAL (z1), 2.5, 1e-9, "z1.real, (2x - 5)^2 = -1");
  188.     gsl_test_rel (GSL_IMAG (z1), 0.5, 1e-9, "z1.imag, (2x - 5)^2 = -1");
  189.   }
  190.  
  191.   {
  192.     gsl_complex z0, z1;
  193.  
  194.     int n = gsl_poly_complex_solve_quadratic (4.0, -20.0, 25.0, &z0, &z1);
  195.  
  196.     gsl_test (n != 2,
  197.           "gsl_poly_complex_solve_quadratic, one root, (2x - 5)^2 = 0");
  198.     gsl_test_rel (GSL_REAL (z0), 2.5, 1e-9, "z0.real, (2x - 5)^2 = 0");
  199.     gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag (2x - 5)^2 = 0");
  200.     gsl_test_rel (GSL_REAL (z1), 2.5, 1e-9, "z1.real, (2x - 5)^2 = 0");
  201.     gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag (2x - 5)^2 = 0");
  202.     gsl_test (GSL_REAL (z0) != GSL_REAL (z1),
  203.           "z0.real == z1.real, (2x - 5)^2 = 0");
  204.     gsl_test (GSL_IMAG (z0) != GSL_IMAG (z1),
  205.           "z0.imag == z1.imag, (2x - 5)^2 = 0");
  206.   }
  207.  
  208.   {
  209.     gsl_complex z0, z1;
  210.  
  211.     int n = gsl_poly_complex_solve_quadratic (4.0, -20.0, 21.0, &z0, &z1);
  212.  
  213.     gsl_test (n != 2,
  214.           "gsl_poly_complex_solve_quadratic, two roots, (2x - 5)^2 = 4");
  215.     gsl_test_rel (GSL_REAL (z0), 1.5, 1e-9, "z0.real, (2x - 5)^2 = 4");
  216.     gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (2x - 5)^2 = 4");
  217.     gsl_test_rel (GSL_REAL (z1), 3.5, 1e-9, "z1.real, (2x - 5)^2 = 4");
  218.     gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (2x - 5)^2 = 4");
  219.   }
  220.  
  221.   {
  222.     gsl_complex z0, z1;
  223.  
  224.     int n = gsl_poly_complex_solve_quadratic (4.0, 7.0, 0.0, &z0, &z1);
  225.  
  226.     gsl_test (n != 2,
  227.           "gsl_poly_complex_solve_quadratic, two roots, x(4x + 7) = 0");
  228.     gsl_test_rel (GSL_REAL (z0), -1.75, 1e-9, "z0.real, x(4x + 7) = 0");
  229.     gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, x(4x + 7) = 0");
  230.     gsl_test_rel (GSL_REAL (z1), 0.0, 1e-9, "z1.real, x(4x + 7) = 0");
  231.     gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, x(4x + 7) = 0");
  232.   }
  233.  
  234.   {
  235.     gsl_complex z0, z1;
  236.  
  237.     int n = gsl_poly_complex_solve_quadratic (5.0, 0.0, -20.0, &z0, &z1);
  238.  
  239.     gsl_test (n != 2,
  240.           "gsl_poly_complex_solve_quadratic, two roots b = 0, 5 x^2 = 20");
  241.     gsl_test_rel (GSL_REAL (z0), -2.0, 1e-9, "z0.real, 5 x^2 = 20");
  242.     gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, 5 x^2 = 20");
  243.     gsl_test_rel (GSL_REAL (z1), 2.0, 1e-9, "z1.real, 5 x^2 = 20");
  244.     gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, 5 x^2 = 20");
  245.   }
  246.  
  247.   {
  248.     gsl_complex z0, z1;
  249.  
  250.     int n = gsl_poly_complex_solve_quadratic (5.0, 0.0, 20.0, &z0, &z1);
  251.  
  252.     gsl_test (n != 2,
  253.           "gsl_poly_complex_solve_quadratic, two roots b = 0, 5 x^2 = -20");
  254.     gsl_test_rel (GSL_REAL (z0), 0.0, 1e-9, "z0.real, 5 x^2 = -20");
  255.     gsl_test_rel (GSL_IMAG (z0), -2.0, 1e-9, "z0.imag, 5 x^2 = -20");
  256.     gsl_test_rel (GSL_REAL (z1), 0.0, 1e-9, "z1.real, 5 x^2 = -20");
  257.     gsl_test_rel (GSL_IMAG (z1), 2.0, 1e-9, "z1.imag, 5 x^2 = -20");
  258.   }
  259.  
  260.   /* Cubic with complex roots */
  261.  
  262.   {
  263.     gsl_complex z0, z1, z2;
  264.  
  265.     int n = gsl_poly_complex_solve_cubic (0.0, 0.0, -27.0, &z0, &z1, &z2);
  266.  
  267.     gsl_test (n != 3, "gsl_poly_complex_solve_cubic, three root, x^3 = 27");
  268.     gsl_test_rel (GSL_REAL (z0), -1.5, 1e-9, "z0.real, x^3 = 27");
  269.     gsl_test_rel (GSL_IMAG (z0), -1.5 * sqrt (3.0), 1e-9,
  270.           "z0.imag, x^3 = 27");
  271.     gsl_test_rel (GSL_REAL (z1), -1.5, 1e-9, "z1.real, x^3 = 27");
  272.     gsl_test_rel (GSL_IMAG (z1), 1.5 * sqrt (3.0), 1e-9, "z1.imag, x^3 = 27");
  273.     gsl_test_rel (GSL_REAL (z2), 3.0, 1e-9, "z2.real, x^3 = 27");
  274.     gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, x^3 = 27");
  275.   }
  276.  
  277.   {
  278.     gsl_complex z0, z1, z2;
  279.  
  280.     int n = gsl_poly_complex_solve_cubic (-1.0, 1.0, 39.0, &z0, &z1, &z2);
  281.  
  282.     gsl_test (n != 3,
  283.           "gsl_poly_complex_solve_cubic, three root, (x+3)(x^2-4x+13) = 0");
  284.     gsl_test_rel (GSL_REAL (z0), -3.0, 1e-9, "z0.real, (x+3)(x^2+1) = 0");
  285.     gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x+3)(x^2+1) = 0");
  286.     gsl_test_rel (GSL_REAL (z1), 2.0, 1e-9, "z1.real, (x+3)(x^2+1) = 0");
  287.     gsl_test_rel (GSL_IMAG (z1), -3.0, 1e-9, "z1.imag, (x+3)(x^2+1) = 0");
  288.     gsl_test_rel (GSL_REAL (z2), 2.0, 1e-9, "z2.real, (x+3)(x^2+1) = 0");
  289.     gsl_test_rel (GSL_IMAG (z2), 3.0, 1e-9, "z2.imag, (x+3)(x^2+1) = 0");
  290.   }
  291.  
  292.   {
  293.     gsl_complex z0, z1, z2;
  294.  
  295.     int n =
  296.       gsl_poly_complex_solve_cubic (-51.0, 867.0, -4913.0, &z0, &z1, &z2);
  297.  
  298.     gsl_test (n != 3,
  299.           "gsl_poly_complex_solve_cubic, three roots, (x-17)^3=0");
  300.     gsl_test_rel (GSL_REAL (z0), 17.0, 1e-9, "z0.real, (x-17)^3=0");
  301.     gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x-17)^3=0");
  302.     gsl_test_rel (GSL_REAL (z1), 17.0, 1e-9, "z1.real, (x-17)^3=0");
  303.     gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (x-17)^3=0");
  304.     gsl_test_rel (GSL_REAL (z2), 17.0, 1e-9, "z2.real, (x-17)^3=0");
  305.     gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, (x-17)^3=0");
  306.   }
  307.  
  308.   {
  309.     gsl_complex z0, z1, z2;
  310.  
  311.     int n =
  312.       gsl_poly_complex_solve_cubic (-57.0, 1071.0, -6647.0, &z0, &z1, &z2);
  313.  
  314.     gsl_test (n != 3,
  315.           "gsl_poly_complex_solve_cubic, three roots, (x-17)(x-17)(x-23)=0");
  316.     gsl_test_rel (GSL_REAL (z0), 17.0, 1e-9, "z0.real, (x-17)(x-17)(x-23)=0");
  317.     gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x-17)(x-17)(x-23)=0");
  318.     gsl_test_rel (GSL_REAL (z1), 17.0, 1e-9, "z1.real, (x-17)(x-17)(x-23)=0");
  319.     gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (x-17)(x-17)(x-23)=0");
  320.     gsl_test_rel (GSL_REAL (z2), 23.0, 1e-9, "z2.real, (x-17)(x-17)(x-23)=0");
  321.     gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, (x-17)(x-17)(x-23)=0");
  322.   }
  323.  
  324.   {
  325.     gsl_complex z0, z1, z2;
  326.  
  327.     int n =
  328.       gsl_poly_complex_solve_cubic (-11.0, -493.0, +6647.0, &z0, &z1, &z2);
  329.  
  330.     gsl_test (n != 3,
  331.           "gsl_poly_complex_solve_cubic, three roots, (x+23)(x-17)(x-17)=0");
  332.     gsl_test_rel (GSL_REAL (z0), -23.0, 1e-9,
  333.           "z0.real, (x+23)(x-17)(x-17)=0");
  334.     gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x+23)(x-17)(x-17)=0");
  335.     gsl_test_rel (GSL_REAL (z1), 17.0, 1e-9, "z1.real, (x+23)(x-17)(x-17)=0");
  336.     gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (x+23)(x-17)(x-17)=0");
  337.     gsl_test_rel (GSL_REAL (z2), 17.0, 1e-9, "z2.real, (x+23)(x-17)(x-17)=0");
  338.     gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, (x+23)(x-17)(x-17)=0");
  339.   }
  340.  
  341.  
  342.   {
  343.     gsl_complex z0, z1, z2;
  344.  
  345.     int n =
  346.       gsl_poly_complex_solve_cubic (-143.0, 5087.0, -50065.0, &z0, &z1, &z2);
  347.  
  348.     gsl_test (n != 3,
  349.           "gsl_poly_complex_solve_cubic, three roots, (x-17)(x-31)(x-95)=0");
  350.     gsl_test_rel (GSL_REAL (z0), 17.0, 1e-9, "z0.real, (x-17)(x-31)(x-95)=0");
  351.     gsl_test_rel (GSL_IMAG (z0), 0.0, 1e-9, "z0.imag, (x-17)(x-31)(x-95)=0");
  352.     gsl_test_rel (GSL_REAL (z1), 31.0, 1e-9, "z1.real, (x-17)(x-31)(x-95)=0");
  353.     gsl_test_rel (GSL_IMAG (z1), 0.0, 1e-9, "z1.imag, (x-17)(x-31)(x-95)=0");
  354.     gsl_test_rel (GSL_REAL (z2), 95.0, 1e-9, "z2.real, (x-17)(x-31)(x-95)=0");
  355.     gsl_test_rel (GSL_IMAG (z2), 0.0, 1e-9, "z2.imag, (x-17)(x-31)(x-95)=0");
  356.   }
  357.  
  358.  
  359.   {
  360.     /* Wilkinson polynomial: y = (x-1)(x-2)(x-3)(x-4)(x-5) */
  361.  
  362.     double a[6] = { -120, 274, -225, 85, -15, 1 };
  363.     double z[6*2];
  364.  
  365.     gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc (6);
  366.  
  367.     int status = gsl_poly_complex_solve (a, 6, w, z);
  368.  
  369.     gsl_poly_complex_workspace_free (w);
  370.  
  371.     gsl_test (status,
  372.           "gsl_poly_complex_solve, 5th-order Wilkinson polynomial");
  373.     gsl_test_rel (z[0], 1.0, 1e-9, "z0.real, 5th-order polynomial");
  374.     gsl_test_rel (z[1], 0.0, 1e-9, "z0.imag, 5th-order polynomial");
  375.     gsl_test_rel (z[2], 2.0, 1e-9, "z1.real, 5th-order polynomial");
  376.     gsl_test_rel (z[3], 0.0, 1e-9, "z1.imag, 5th-order polynomial");
  377.     gsl_test_rel (z[4], 3.0, 1e-9, "z2.real, 5th-order polynomial");
  378.     gsl_test_rel (z[5], 0.0, 1e-9, "z2.imag, 5th-order polynomial");
  379.     gsl_test_rel (z[6], 4.0, 1e-9, "z3.real, 5th-order polynomial");
  380.     gsl_test_rel (z[7], 0.0, 1e-9, "z3.imag, 5th-order polynomial");
  381.     gsl_test_rel (z[8], 5.0, 1e-9, "z4.real, 5th-order polynomial");
  382.     gsl_test_rel (z[9], 0.0, 1e-9, "z4.imag, 5th-order polynomial");
  383.   }
  384.  
  385.   {
  386.     /* : 8-th order polynomial y = x^8 + x^4 + 1 */
  387.  
  388.     double a[9] = { 1, 0, 0, 0, 1, 0, 0, 0, 1 };
  389.     double z[8*2];
  390.  
  391.     double C = 0.5;
  392.     double S = sqrt (3.0) / 2.0;
  393.  
  394.     gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc (9);
  395.  
  396.     int status = gsl_poly_complex_solve (a, 9, w, z);
  397.  
  398.     gsl_poly_complex_workspace_free (w);
  399.  
  400.     gsl_test (status, "gsl_poly_complex_solve, 8th-order polynomial");
  401.  
  402.     gsl_test_rel (z[0], -S, 1e-9, "z0.real, 8th-order polynomial");
  403.     gsl_test_rel (z[1], C, 1e-9, "z0.imag, 8th-order polynomial");
  404.     gsl_test_rel (z[2], -S, 1e-9, "z1.real, 8th-order polynomial");
  405.     gsl_test_rel (z[3], -C, 1e-9, "z1.imag, 8th-order polynomial");
  406.     gsl_test_rel (z[4], -C, 1e-9, "z2.real, 8th-order polynomial");
  407.     gsl_test_rel (z[5], S, 1e-9, "z2.imag, 8th-order polynomial");
  408.     gsl_test_rel (z[6], -C, 1e-9, "z3.real, 8th-order polynomial");
  409.     gsl_test_rel (z[7], -S, 1e-9, "z3.imag, 8th-order polynomial");
  410.     gsl_test_rel (z[8], C, 1e-9, "z4.real, 8th-order polynomial");
  411.     gsl_test_rel (z[9], S, 1e-9, "z4.imag, 8th-order polynomial");
  412.     gsl_test_rel (z[10], C, 1e-9, "z5.real, 8th-order polynomial");
  413.     gsl_test_rel (z[11], -S, 1e-9, "z5.imag, 8th-order polynomial");
  414.     gsl_test_rel (z[12], S, 1e-9, "z6.real, 8th-order polynomial");
  415.     gsl_test_rel (z[13], C, 1e-9, "z6.imag, 8th-order polynomial");
  416.     gsl_test_rel (z[14], S, 1e-9, "z7.real, 8th-order polynomial");
  417.     gsl_test_rel (z[15], -C, 1e-9, "z7.imag, 8th-order polynomial");
  418.  
  419.   }
  420.  
  421.   {
  422.     int i;
  423.  
  424.     double xa[7] = {0.16, 0.97, 1.94, 2.74, 3.58, 3.73, 4.70 };
  425.     double ya[7] = {0.73, 1.11, 1.49, 1.84, 2.30, 2.41, 3.07 };
  426.  
  427.     double dd_expected[7] = {  7.30000000000000e-01,
  428.                                4.69135802469136e-01,
  429.                               -4.34737219941284e-02,
  430.                                2.68681098870099e-02,
  431.                               -3.22937056934996e-03,
  432.                                6.12763259971375e-03,
  433.                               -6.45402453527083e-03 };
  434.  
  435.     double dd[7], coeff[7], work[7];
  436.     
  437.     gsl_poly_dd_init (dd, xa, ya, 7);
  438.  
  439.     for (i = 0; i < 7; i++)
  440.       {
  441.         gsl_test_rel (dd[i], dd_expected[i], 1e-10, "divided difference dd[%d]", i);
  442.       }
  443.  
  444.     for (i = 0; i < 7; i++)
  445.       {
  446.         double y = gsl_poly_dd_eval(dd, xa, 7, xa[i]);
  447.         gsl_test_rel (y, ya[i], 1e-10, "divided difference y[%d]", i);
  448.       }
  449.  
  450.     gsl_poly_dd_taylor (coeff, 1.5, dd, xa, 7, work);
  451.     
  452.     for (i = 0; i < 7; i++)
  453.       {
  454.         double y = gsl_poly_eval(coeff, 7, xa[i] - 1.5);
  455.         gsl_test_rel (y, ya[i], 1e-10, "taylor expansion about 1.5 y[%d]", i);
  456.       }
  457.   }
  458.  
  459.  
  460.   /* now summarize the results */
  461.  
  462.   exit (gsl_test_summary ());
  463. }
  464.